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ABSTRACT 

We present the first 3D Monte Carlo (MC) photoionisation code to include a fully self- 
consistent treatment of dust radiative transfer (RT) within a photoionised region. This is the 
latest development (Version 2.0) of the gas-only photoionisation code MOCASSIN (Ercolano 
et al., 2003a), and employs a stochastic approach to the transport of radiation, allowing both 
the primary and secondary components of the radiation field to be treated self-consistently, 
whilst accounting for the scattering of radiation by dust grains mixed with the gas, as well as 
the absorption and emission of radiation by both the gas and the dust components. An escape 
probability method is implemented for the transfer of resonance lines that may be absorbed 
by the grains, thus contributing to their energy balance. The energetics of the co-existing dust 
and gas components must also take into account the effects of dust-gas collisions and pho- 
toelectric emission from the dust grains, which are dependent on the grain charge. These are 
included in our code using the average grain potential approximation scheme. 

A set of rigorous benchmark tests have been carried out for dust-only spherically symmet- 
ric geometries and 2D disk configurations, mocassin's results are found to be in agreement 
with those obtained by well established dust-only RT codes that employ various approaches 
to the solution of the RT problem. 

A model of the dust and of the photoionised gas components of the planetary nebula (PN) 
NGC 3918 is also presented as a means of testing the correct functioning of the RT procedures 
in a case where both gas and dust opacities are present. The two components are coupled via 
the heating of dust grains by the absorption of both UV continuum photons and resonance line 
photons emitted by the gas. The MOCASSIN results show agreement with those of a ID dust 
and gas model of this nebula published previously, showing the reliability of the new code, 
which can be applied to a variety of astrophysical environments. 
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1 INTRODUCTION 



The presence of dust grains in ionised plasma environments can 
have significant effects on the radiative transfer (RT) and influence 
the physical conditions of the gas, as the grains compete with the 
gas for the absorption of continuum UV photons. As well as be- 
ing heated via the absorption of such photons, the dust grains are 
also heated by nebular resonance line radiation emitted by the gas. 
Moreover, photo-electric emission from dust grains and dust-grain 
collisions provide additional cooling and heating channels for both 
the gas and dust components. The evident coupling between the 
co-existing dust and gas phases can only be treated properly in a 
photoionisation code by the incorporation, in the RT, of scattering, 



absorption and emission of radiation by dust particles mixed with 
the gas. 

Whilst some ID photoionisation codes, such as CLOUDY (Fer- 
land et al. 1998, Van Hoof et al. 2004) and the Harrington code 
(Harrington et al., 1988, hereafter H88), include a treatment for 
the RT of dust, and a number of 2D or 3D pure-dust RT codes us- 
ing MC or ray tracing techniques are available (e.g. Bjorkman & 
Wood, 2001, or Pascucci et al., 2004 for a recent review), to our 
knowledge, no existing 3D photoionisation code self-consistently 
treats the transfer of both the stellar and diffuse components of the 
radiation field through a region where dust particles are mixed with 
the gas. 

We present here a new version (2.0) of the 3D photoionisa- 
tion code MOCASSIN (Ercolano et al., 2003a, hereafter Paper I), 
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designed to achieve the above by employing a stochastic approach 
to simultaneously solve the dust and gas RT. As well as calculating 
emission line spectra and the 3D ionisation and electron temper- 
ature and density structure of a nebula, the code can also deter- 
mine accurate dust temperatures and spectral energy distributions 
(SEDs) by treating discrete grain sizes and different grain species 
separately, compared to the frequently used approximation of hav- 
ing a single grain as representative of an ensemble of grains. 

The new code is described here together with its solutions 
for a number of pure-dust benchmark problems and the pseudo- 
benchmark case of a dust and gas model for the PN NGC 3918. 
Section 2 contains a description of the MC techniques applied to 
the solution of the RT problem in a dusty photoionised region. 
Dust-only benchmark models and their solutions are presented in 
Section 3 for the ID and the 2D cases. A 3D model of the pho- 
toionised region and the thermal IR emission of the PN NGC 3918 
is presented in Section 4, together with a comparison of our results 
with those obtained by previous ID modelling and the available 
observational data. A final summary is presented in Section 5. 



2 DUST RADIATIVE TRANSFER IN THE 
PHOTOIONISED REGION 

Paper I presented the fully 3D photoionisation code, MOCASSIN, 
which employs a MC approach to the transfer of radiation. The de- 
scription of the radiation field as composed of a discrete number of 
monochromatic packets of energy Eo (Abbot & Lucy, 1985) allows 
for both primary and secondary components of the radiation field 
to be determined exactly. The mean intensity of the radiation field, 
J„, is derived at each location using the MC estimator constructed 
by Lucy (1999) as given by equation 1 1 of Paper I. The path of an 
energy packet through the grid completely defines its contribution 
to J„ at each location; this reduces the RT problem to the compu- 
tation of the energy packets' trajectories as they diffuse out of the 
nebula. In the case of a pure-gas nebula, only absorption and re- 
emission of the packets by the ions in the gas need be considered, 
and therefore the local gas opacities and emissivities, which depend 
on the electron temperatures and ionic fractions, are the necessary 
ingredients of the computation. 

However, dust grains that are mixed with the gas in the ionised 
region provide an extra source of opacity and, as well as competing 
for the absorption of UV continuum radiation, are also capable of 
absorbing line photons emitted by resonant transitions. In the re- 
mainder of this section the discussion will be limited to the transfer 
of the continuum radiation, with the description of the resonant line 
transfer escape probability method being given in Section l2.3.1l 

The random walk of stellar energy packets injected into the 
nebula is characterised by absorption and re-emission events due 
to the atoms and ions in the gas, as well as absorption, re-emission 
and scattering events due to the dust grains. The locations of the in- 
teractions are calculated using the (second) MC method described 
in Section 2.4 of Paper I, with the random path of a packet between 
events, I, being described by 
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where t„ p is the total optical depth of the gas and dust to the next 
interaction at frequency vp along the packet's direction of travel, 
/vgp" is the local gas opacity at up, while Kl xi = Kj ca +Kj bs is the 
dust extinction opacity at up due to scattering and absorption by 
the dust grains. 

This first MC choice can only tell us whether a packet is going 



to interact at a given location or is going to continue into the next 
grid cell. The nature of the interaction is however not known at this 
stage. The probability of a packet of frequency up undergoing a 
scattering event, P SC a(up), is given by the ratio of the scattering 
to the total opacity, comprising both the gas and dust extinction 
opacities: 
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(2) 



If the second MC test, performed at this stage, results in the packet 
being absorbed, then the process continues as for the gas-only ver- 
sion of the code described in Paper I. 

Energy conservation is enforced by ensuring that the packets 
that are absorbed by the dust grain are re-emitted in the new di- 
rection immediately after the absorption event. The frequency of 
the re-emitted packet is selected according to the local gas+dust 
emissivities, with the frequency distribution of the re-emitted pack- 
ets determined from the normalised cumulative probability density 
function. A grain of size a and species s will provide a contribution 
to the local continuous emission at frequency u equal to 

ji7\v)=^aQ't:{u)B v {T^) (3) 

where Qa b s s (V) is the grain absorption efficiency at frequency u and 
Bi/(7a,s) is the Planck function at the grain temperature. If j sas {u) 
is the frequency-dependent gas continuum emissivity, corrected for 
the contribution due to He I and He II Lyman lines, as described 
in Section 2.6 of Paper I, then the total emission from a volume 
element at frequency u is given by 
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(4) 



where A r gas is the gas number density, n ajS is the number density 
of grains of species s and size a, while the summation and the in- 
tegration in the second term of the RHS are over the grain species 
included and the size distribution, respectively. 

The normalised cumulative probability density function, p„, 
is then calculated as follows 



p v = 



(5) 



where the integral in the denominator is over the full spectral range 
considered and defined by u m i n and f mas . 

A packet undergoing a scattering event will simply be re- 
directed, with its frequency remaining unchanged. Unless a phase 
function is specified, scattering is assumed to be isotropic, with 
the new random directions being selected according to equation 4 
in Paper I. In the non-isotropic case the Henyey-Greenstein (HG) 
phase function (Henyey & Greenstein, 1941) can be used in its 
standard or modified form (Henney & Axon, 1995) to determine 
a packet's scattering angle. In this case the anisotropy parameter g 
must be specified by the user. 

A packet is transferred through the grid until no more interac- 
tions with gas or dust are possible. In general, emission lines are 
assumed to escape without further interaction; special attention is 
given to resonance lines, which can be absorbed by the dust grains, 
as we discuss later in Sectionl2.3.1l 



2.1 The dust opacities 

MOCASSIN uses standard Mie scattering theory (e.g. Kattawar & 
Plass 1967; Wickramasinghe 1972; Andriesse 1979) to evaluate 
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the effective absorption and scattering efficiencies, Q a bs{o., A) and 
Qsca{a, A), for a grain of radius a at wavelength A, from the di- 
electric constants of any given material. The optical constants li- 
brary distributed with the source code includes the most frequently 
used astronomical dust species, such as warm and cold silicates 
(Ossenkopf, Henning & Mathis 1992), astronomical silicates and 
graphite (Draine & Laor 1993, Draine & Lee 1984, Hageman, 
Gudat & Kunz 1974), amorphous carbon (Hanner 1988), aSiC 
(Pegourie 1988) and another seventeen species. However, supple- 
mentary/alternative data can be easily supplied by the user. 

With the exception of silicates, graphite and SiC by Draine 
& Laor (1993) and graphite by Hageman, Gudat & Kunz (1974), 
the data provided does not extend to the far UV (FUV) domain. At 
ionising energies, such data must be supplemented in order to allow 
for the absorption of starlight that is the main heating process for 
the grains, which in turn may provide a significant source of opacity 
to the nebula. 



and neglecting photoelectric emission processes and gas-dust colli- 
sions (see Section l2.3.2t . the radiative equilibrium equation can be 
written as follows: 

f B x (T)Q abs (a,X)d\ = (7) 
Jx(T) Q abs (a, \)d\ + ^2H L 

-min L 

where B\(T) is the usual Planck function evaluated at the grain 
temperature and integrated over the frequency range delimited by 
v m in and v m ax and the summation in the second term of the right- 
hand-side of the equation is over all the resonance lines included in 
the model. 



2.2 Grain sizes and Chemistry 

Grains of different sizes and species are treated separately, with in- 
dividual grain temperatures being determined for each grain radius 
and species included. It has been previously demonstrated that ap- 
proximating a grain mixture with a representative composite grain 
can yield misleading conclusions (e.g. Van Hoof et al. 2004), given 
that the emission due to the smaller, hotter grains may be hard to 
reproduce accurately using this technique. 

There is no limit to the number of sizes that may be included 
in a given model, except perhaps that intrinsically imposed by the 
available computing resources. Any type of grain size distribution 
and chemical mixture can be handled by the code. Furthermore, 
these do not need to be defined homogeneously across the grid, but 
multiple chemistry and/or grain-size distributions can be specified 
at different locations. The latter option was implemented in view 
of the dual circumstellar dust chemistry (oxygen- and carbon-rich) 
found in some PNe ionised by Wolf-Rayet central stars (e.g. Cohen 
et al., 2002; De Marco et al. 2004). 

2.3 Determination of dust temperatures 

The local grain temperatures for each species and grain size are de- 
termined by the equation of radiative equilibrium between all cool- 
ing and heating channels. The grains are mainly cooled via their 
temperature-dependent thermal emission and heated by the absorp- 
tion of UV photons from the radiation field, as well as by reso- 
nance line photons. As will be described in more detail in the next 
section, resonance line transfer is treated using an escape probabil- 
ity method similar to that described by Cohen, Harrington & Hess 
(1984) and later modified and applied by H88 to their dust model 
for the PN NGC 3918. At each location in the grid, we estimate 
/ c , the position-dependent fraction of photons in a given resonance 
line that escape, and, if Gl is the rate of energy generation in line 
L, then, for any given dust species, the heating contribution from 
line L to a grain of radius a is given by 

H L =G L -(1- fe) ■ Qabs{a, A L )/ Kd (A L ) (6) 

where Al is the central wavelength of the resonance line and 
td(Ai,) is the mean dust opacity of that species, which determines 
the intensity within the line. 

The continuous part of the mean intensity of the radiation 
field, J(A), is estimated as discussed in Section|2j using the above 



2.3.1 Resonant line transfer 

In the previous gas-only version of the MOCASSIN code, line pack- 
ets could escape immediately after they had been created. The only 
exception was made in the case of He I Lyman lines, that can fur- 
ther ionise H°, and He II Lyman lines that can ionise both H° and 
He . In the presence of dust grains, however, care must also be 
taken in the transfer of resonance emission lines such as, for exam- 
ple, H I Lya, C IV A1549, N v A1240, C II A1336, Si IV A1397 
and Mg II A2800. These lines, amongst others, can undergo many 
resonant scattering events, resulting in a much longer random walk 
before they can reach the edge of the nebula and escape, with a 
consequently higher probability of being absorbed by dust. 

The emergent resonant line intensities will, therefore, be atten- 
uated by dust absorption, with the absorbed fraction contributing to 
the heating of the grains. This is an extra complication arising from 
the interaction between the gas and the dust phase, it is nevertheless 
essential to carefully treat this process, which can often have siz- 
able effects on the emergent continuous spectral energy distribution 
by raising the grain temperatures, as well as affecting the observed 
intensities of the emergent resonance emission lines. 

H88 employed a generalisation to ID spherical geometry 
of the plane-parallel escape probability approach of Cohen et al. 
(1984), whereby the resonance line photons are assumed either to 
be absorbed locally by the dust or be scattered into the line wings 
and escape in a single flight. We have further generalised the argu- 
ment to arbitrary geometry and implemented it within the 3D MO- 
CASSIN code, for the transfer of some important resonance emis- 
sion lines, including Lya 1216, C II 1336, C IV 1549, N V 1240, 
Mg II 2800 and Si IV 1397. Any other line of interest can easily be 
added to this list by providing the relevant atomic data to the code's 
library. 

The basic idea consists of obtaining the fraction of escaping 
resonance line photons, / c , at each grid location for each resonance 
line. This quantity can then be used in equation|S|above in order to 
estimate the contribution of a given resonance line to the dust heat- 
ing. The contribution from this grid location to the dust-attenuated 
nebular luminosity in this resonance line can subsequently be ob- 
tained by multiplying this factor by the value calculated using the 
formal solution or the full MC method, as described in Paper I. 

For each resonance line P c [t(u)], the probability of escape 
without further interaction of a photon scattered in direction u, is 
first calculated. Values of P e [r(u)] are calculated at each location 
in the grid, in a number of radial directions emerging from the cen- 
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tre of the cell 1 using the following expression (Kwan & Krolik, 
1981): 



-0F(1.2 + b) 



P e [r(u)] = 



2t 



T > 1 



r<l 



(8) 



where fe = yJlog(r)/ (1 + t/t w ), r = t(«) is the optical depth at 
line centre in direction u and t w = 10 5 . Assuming the scattering to 
be isotropic, we obtain P c = < P c [t(u)] >, the escape probability 
averaged over all directions. Using this value of P c and equation 
A 8 of Cohen et al. (1984) we can estimate the total fraction of es- 
caping resonance line photons, / e . It is perhaps worth noting at this 
point that our approach differs somewhat from that of H88, who 
used equation Al of Cohen et al. (1984) to calculate Mi [t(w)], the 
dust-free escape probability for a resonance line photon travelling 
in direction u, and then obtained a direction averaged escape prob- 
ability P e , by taking the mean of the M2 [t(m)] values for 40 radial 
directions. 

It should, finally, be stated that in a dust-only simulation the 
only source of grain heating comes from the continuous radiation 
field, so that in this case the above becomes irrelevant and no over- 
heads are introduced in the execution time/memory requirements. 

We compared the treatment described above with the exact 
results of Hummer & Kunasz (1980) by constructing a plane paral- 
lel model in MOCASSIN to mimic a uniform slab of optical thick- 
ness T. TableQshows the results obtained by MOCASSIN and those 
given in Table 2 of Hummer & Kunasz (1980) for the fraction of 
escaping resonance line photons from mid-plane. We included the 
cases whenT = 10 4 and 10 6 , and for a range of f3 values, where /3 
is the ratio of continuum to line opacity, Kd/ni. For the T = 10 4 
case we also include the results given in Table 4 of Cohen et al. 
(1984). All our values are within 20% of the exact values, apart 
from the T = 10 6 , f3 = 4.03 x 10~ 7 case, where a discrepancy of 
34% is obtained. The results of Hummer & Kunasz (1980) depend 
also on the parameter a, which is the ratio of natural to doppler 
line widths. Only the latter is taken into account by our treatment 
and that of Cohen et al. (1984), therefore our results are best com- 
pared with those listed by Hummer & Kunasz (1980) for the lowest 
value of a (4.7 x 10 -4 ). We note that although we have reasonable 
agreement with the results of Hummer & Kunasz (1980) for the 
a = 4.7 x 10~ 4 and T < 10 6 case, the neglect of the damping 
wings for transitions with large optical depth means that we will 
necessarily overestimate the fraction of photons escaping for larger 
values of a and T. 



2.3.2 Additional heating and cooling channels 

A number of additional microphysical processes come into play 
when dust and gas interactions are considered. Photoelectric emis- 
sion from the surface of dust grains can provide an additional heat- 
ing channel for the photoionised gas in a variety of astrophysical 
environments, including H II regions and PNe (e.g. Maciel & Pot- 
tasch 1982, Baldwin et al. 1991, Borkowski & Harrington 1991, 
Ercolano et al., 2003c). When this process is included, the energy 
carried out by the photoelectrons must also be accounted for in the 
thermal balance equations of the grains as a heat sink. The calcula- 
tion of photoelectric thresholds and yields requires a knowledge of 



64 directions were used at each cell for the model presented in SectionEI 



Table 1. Resonance line photon escaping fractions from a slab geometry of 
optical thickness T= 10 4 and 10 6 , calculated for a range of continuum- to 
line-opacity ratios, /? = kJki . MOCASSIN'S results for the mid-plane values 
(MOC) are compared to those obtained by Cohen et al. (1984) (CHH84, for 
T= 10 4 ) and to the exact treatment of Hummer & Kunasz, 1980 (HK80). 
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the individual grain potentials (Weingartner & Draine 2001). Fur- 
thermore, energy is also exchanged between the gas and the grains 
via collisions that will cool the gas and heat the grains in the case 
when the latter have the lower temperature. 

The processes listed above are treated by MOCASSIN using 
the average grain potential approximation, as described by Baldwin 
et al. (1991), that is based on the detailed balance between grain 
charge gain and loss rates. Whilst this is an excellent approximation 
for large grains, it may fail to provide an accurate description for 
very small grains (e.g. Weingartner and Draine, 2001). 

As photoelectric emission from dust grains and gas-grain col- 
lisional processes are irrelevant to the applications presented in this 
work, we refer to a forthcoming paper (Ercolano et al., in prepara- 
tion) for a more detailed description of their implementation in the 
MOCASSIN code. 

Finally, we note that at present grain charge transfer is not 
included in our code. This may provide a net recombination process 
for the heavy elements (e.g. Draine & Sutin 1987) in the H + region. 



2.4 Emerging Spectral Energy Distribution 

The emergent SED is obtained by collecting the energy packets as 
they reach the outer edge of the grid and escape. Frequency and 
direction information is retained for each packet, so that it is later 
possible to account for orientation effects, that are often of vital 
importance for asymmetric models. 

As mentioned above, the spectral distribution of the radiation 
re-processed by the grains is calculated by integrating equation 3 
over the size distribution and abundances of the grain mixture, with 
B v (T a ,s) being evaluated for each grain of species s and radius 
a. This means that the emerging continuum energy packets will 
automatically provide the correct spectral shape for the chosen dust 
chemistry and discrete size distribution chosen in the model. 

This is certainly an improvement over the practice of calcu- 
lating SEDs using a single dust temperature to represent an ensem- 
ble of grains of various species and sizes. Whilst this may have 
been a necessity in the past, due to limited computing power, it 
may result in misleading conclusions, as the emission from hotter, 
smaller grains may not be accurately reproduced (e.g. Van Hoof et 
al. 2004). 

Furthermore, very small grains that have a small heat capacity 
can undergo large temperature fluctuations, also known as temper- 
ature spikes. This process, first suggested by Greenberg (1968), is 
included in the current version of MOCASSIN, using the method of 
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Table 2. Input parameters for spherically symmetric benchmark models. 
For all models, the ratio of the outer radius i? ou t to the inner radius of the 
shell Y = 1000; N dust is the number density of the grain; L t and T c ff are 
the luminosity and the effective temperature of the ionising source, which 
radiates as a blackbody; ao is the radius of the dust grains. 
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Guhathakurta & Draine (1989) in the implementation of Sylvester 
et al. (1997). This effect is generally important for grains with typ- 
ical radii of approximately 10 A, however the minimum grain sizes 
for which the process should be taken into account is treated as an 
input parameter in order to provide the user with some control over 
the overheads introduced. 



3 PURE DUST BENCHMARKS 

mocassin's dust and gas processes are fully integrated, it is how- 
ever possible to run the code for models containing only one or the 
other types of process. The pure gas photoionisation version was 
benchmarked against established ID photoionisation codes accord- 
ing to the Lexington-Meudon test models (Paper I), and we have 
checked that, in the absence of dust, these can still be correctly 
reproduced by the current version. Furthermore, the new dust RT 
procedures were thoroughly tested by running dust-only simula- 
tions for a set of ID and 2D benchmark problems, as described 
below. 

As for previous versions, MOCASSIN Version 2.0 is also fully 
parallel and can run efficiently on multi-processors or Beowulf sys- 
tems. Dust-only models, however, converge much more quickly 
than full photionisation ones, thanks the to dust opacities being 
virtually independent of temperature. The amount of memory re- 
quired to run such models is also much less, which makes MO- 
CASSIN also suitable for running on smaller single-processor ma- 
chines. As a guideline, all benchmark models presented in this sec- 
tion were computed on a Pentium M 1.70 GHz portable PC with 
1GB of RAM. Models converged to acceptable accuracy within 
1 hr or less, with longer computational times sometimes being re- 
quired to smooth out numerical errors that are intrinsic to the MC 
technique employed. 

3.1 Spherically symmetric benchmarks 

A set of spherically symmetric benchmarks problems have been 
provided by Ivezic et al. (1997), where solutions obtained by three 
RT codes, including the widely used DUSTY (Ivezic & Elitzur, 
1997), are compared. The codes, each implementing a different nu- 
merical scheme, were found to agree to better than 0. 1 per cent. The 
radial dust temperature distributions and the emerging SEDs for a 
number of these benchmark tests were calculated using MOCASSIN 
and compared to results obtained with the current version of DUSTY 
(2.01). Some preliminary benchmark results have already been pre- 
sented by Ercolano et al. (2005). 

We present here three spherically symmetric homogeneous 
dust shell models, illuminated by a central blackbody source and 



t (1 M m)=1, 10, 100 




Figure 1. Spherically symmetric, homogeneous benchmark models for a 
10 4 Lq star surrounded by dust shells with Ti=l, 10 and 100. The thicker 
lines have been used to represent the more optically thick models. The SEDs 
predicted by DUSTY (dashed lines) are compared to those obtained by MO- 
CASSIN (solid line). 



Table 3. Input parameters for 2D disk benchmark models, tv: mid-plane 
visual optical depth (A = 55()nm); i?; n : inner disk radius; Rout - outer disk 
radius; z d : disk thickness; L„[Lq]: central star luminosity; T c g [K]: central 
star effective temperature; ao[/im]: grain radius. The density distribution 
was calculated using the parameters listed and equation 4 of Pascucci et al. 
(2004). 
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containing silicate grains of radius 0. 16/im. Optical constants for 
the 'astronomical silicates' of Draine & Lee (1984) were used in 
the MOCASSIN models, as these are also included in DUSTY's stan- 
dard optical constants library. The scattering was assumed to be 
isotropic for all benchmark models. The optical thickness of the 
shell at A = 1 /im for the three models was n=l, 10 and 100. The 
input parameters are summarised in Table|2| 

FigureQshows a comparison between the SEDs predicted by 
MOCASSIN (solid lines) and those predicted by DUSTY (Ivezic et 
al., 1997; dashed line) for the n = 1 (thin lines), n = 10 (medium 
line) and n = 100 (thick line) models. The emerging SEDs are 
plotted, by analogy with Figure 2 of Ivezic et al. (1997), as di- 
mensionless, distance- and luminosity-independent spectral shapes 
\F\I J F\d\. Taking into account that a completely different ap- 
proach to the radiative transfer and the calculation of the SEDs is 
used by the two codes, the results obtained from the benchmark- 
ing are in very good agreement. The radial distributions of grain 
temperatures (not shown here) were also found to be in excellent 
agreement. 
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Figure 2. Spectral energy distributions for disk inclinations i = 12.5° (left panel) and i = 77.5° (right panel) for the Ty = 0. 1 (solid line), ry = 1 (dotted line), 
tv = 10 (dashed line) and ry = 100 (dot-dash line) benchmark models around a 5800 K, 1 L0 . This figure is directly comparable to figure 7 of Pascucci et 
al. (2004). 



3.2 2D disk benchmarks 

Four benchmark solutions for the continuum RT in a 2D disk con- 
figuration were presented by Pascucci et al. (2004). They tested 
the behaviour of five different codes for a set of well-defined 2D 
models having different visual (A=550nm) optical depths, tv, in- 
cluding a high optical depth (tv=100) and strong (isotropic) scat- 
tering case. All models consisted of a point-source embedded in a 
circumstellar dust disk with an inner dust-free cavity. The dust is 
composed of spherical astronomical silicate grains (Draine & Lee, 
1984) with a radius of 0.12/^m and density of 3.6gcm™ 3 . Results 
from various codes for the radius-dependent grain temperatures and 
the emergent SEDs for a number of radial directions and disk in- 
clinations were compared by Pascucci et al. (2004) and found to 
be in reasonable agreement for most test cases. All benchmark 
tests were also successfully reproduced by the MOCASSIN code as 
briefly shown in this section. The disk geometry employed is typ- 
ical of those generally used for the study of disks around T Tauri 
stars (Natta et al. 2000) and described by equation 4 of Pascucci et 
al. (2004). A summary of all input parameters is given in Table|3| 
with the symbols used being defined in the caption. 

Figure l3~2l shows the SEDs (A F\ [W m~ 2 ] against A [/im]) for 
two disk inclinations: one almost face-on (i = 12.5° - left panel), 
and the other almost edge-on (i = 77.5° - right panel -) for the 
tv = 0.1 (solid line), 1 (dotted line), 10 (dashed-line) and 100 (dash- 
dot line) benchmark models around a 5800 K 1 L0 star. This figure 
is directly comparable with figure 7 of Pascucci et al. 2004 and 
shows good agreement for all cases. 

The published benchmark solutions for the radial dust temper- 
ature distributions were also successfully matched by MOCASSIN. 



4 MODELLING THE THERMAL EMISSION OF THE PN 
NGC 3918 

The benchmark models described above were designed to test the 
solution of the dust RT problem in the absence of gas. A set of 
benchmarks to test the correct RT solution of dust and gas mixed 
together is however not available. In this work we will present 
a fully-self consistent 3D photoionisation model of the PN NGC 
3918, where dust grains are mixed with the gas inside the ionised 



region. The geometry of this model and the detailed dust proper- 
ties are based on the work of Clegg et al. (1987) and H88, who 
used a volume weighted combination of two spherically symmet- 
ric dust and gas density distributions fed into a ID photoionisation 
code that could also account for dust absorption. Their results are 
treated here as a 'realistic' benchmark case that can be compared 
with MOCASSIN'S solution, within the limits imposed by the differ- 
ent approaches employed by the two codes and the different atomic 
data sets available to them. 



4.1 Photoionisation models of NGC 3918 

A pure-gas photoionisation model of NGC 39 1 8 was originally pre- 
sented by Clegg et al. (1987), where the nebula was approximated 
by two dense (optically thick) cones immersed in a sphere of more 
diffuse (optically thin) gas. This biconical model, which was suc- 
cessful in reproducing most of the spectroscopic constraints avail- 
able for this object, became unrealistic once the morphology of 
NGC 3918 had been disclosed by later observations (e.g. Corradi 
et al., 1999). Nevertheless, the same density distribution was used 
by Ercolano et al. (2003b) to construct a 3D model of NGC 3918 
using MOCASSIN; the model was used as a 'realistic' 3D bench- 
mark for the then new photoionisation code. Furthermore, differ- 
ences between the two sets of results for the ionisation structure and 
emission line spectra were identified in that paper and discussed in 
the light of the effects of the radiation field being transferred self- 
consistently through the discontinuity between the cones and the 
surrounding material. A new model for the pure-gas photoionised 
region of NGC 3918 was also presented by Ercolano et al (2003b). 
This provided a better fit to the available spatio-kinematical con- 
straints, as well as reproducing most of the observed spectroscopic 
characteristics. However, in the present work, NGC 3918 is treated 
merely as a benchmark model to test the coupling of dust and gas 
in MOCASSIN'S radiative transfer and, for this purpose, only the 
original bi-conical model needs to be considered. 



4.2 The dust model 

The thermal IR emission of NGC 3918 was modelled by H88, who 
added dust grains to the ID photoionisation model described above. 
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Table 4. Input parameters for the MOCASSIN dust and gas model of 
NGC 3918. All gas abundances are given by number. See text for a de- 
scription of the symbols used. 



L» 


7224 L Q 


N/H 


1.5x10" 


-4 


T e ff 


140000 K 


O/H 


5.0x10" 


-4 


P 


2.2 gem -3 


Ne/H 


1.2x10" 


-4 


&min 


0.04^m 


Mg/H 


1.4x10" 


-5 


d-max 


0.40^m 


Si/H 


1.0x10" 


-5 


R 


0.0005 


S/H 


1.6x10" 


-5 


He/H 


0.107 


Ar/H 


2.0x10- 


-6 


C/H 


8.0xl0~ 4 


Fe/H 


3.7x10" 


-7 



They explored a number of size distributions and dust chemistries 
and agreed on a best-fitting model, which they named the 'pri- 
mary model', consisting of spherical graphite grains, with p — 
2.2 gem -3 (Draine & Lee, 1984 and Hagemann, Gudat & Kunz, 
1974) with the standard MRN dust parameters of a m i n =0.04fj,m 
and a ma x=0.40/im and a dust-to-hydrogen ratio, R, of 0.0005 by 
mass. 

The same input parameters as for the pure-gas photoionisation 
model were used for the gas and dust model, except for the lumi- 
nosity of the central star, which was increased by 4.7% in order 
to compensate for the stellar radiation absorbed by the dust (H88 
argued that the same result could have also been obtained by de- 
creasing the gas density by 2.3%). Due to the total optical depth of 
the dust being less than unity in this model, the effects of the inclu- 
sion of the dust opacity on the emission line spectrum, ionisation 
and gas thermal structure of the nebula are small. Results for all 
the gas properties were found to be similar to those obtained in the 
dust-free models, as described by Ercolano et al. (2003b). 

The central star parameters and the elemental abundances used 
for the dust and gas 3D MOCASSIN model are summarised in Ta- 
ble|4] the parameters controlling the shape of bicone were given in 
Table 1 of Ercolano et al. (2003b). 

It is worth noting at this point that the effects of the scattering 
opacity on the radiation field were not taken into account by H88. 
For ease of comparison, this process was also neglected in the ini- 
tial MOCASSIN models, to be reactivated in the subsequent models, 
in order to investigate its effects on the final results. We found that 
for the relatively low dust optical depths of NGC 3918, dust scat- 
tering had only a small effect on the final dust temperatures and 
SEDs. 

Photoelectric emission from dust grains and gas-grain colli- 
sions were not included in the models of H88 and, therefore, were 
also neglected in our models. Finally, as the MRN size distribu- 
tion was truncated to a lower limit of 0.04/im, the grains are not 
expected to be affected by temperature spiking (see section 2.4). 

4.3 Comparison with previous models and available IR 
observations 

The observational IR data available for NGC 3918, as described 
by H88, include four sets in the wavelength range 8 to 100 /im. 
Broadband photometry at 12, 25, 60 and 100/im is available from 
the IRAS Point Source Catalog (1985), and at 36 and 70 /im from 
Moseley (1980). Cohen & Barlow (1980) gave narrow band pho- 
tometry from 8/im to 20^im, while a similar spectral range (8/im to 
23/im) was also covered with the IRAS Low-Resolution Spectrom- 
eter (Pottasch et al., 1986). We refer to H88 and to the respective 
references listed above for further description of these data sets. 



Table 5. IR fluxes [Jy] predicted by the MOCASSIN model and by H88, 
compared with line- and coloured-corrected IRAS Point Source Catalog 
fluxes and the 9.6/xm flux measurement of Cohen & Barlow (1980) 





9.6fim 


12^tm 


25 (im 


60/im 


100/im 


Observed 


1.04 


4.7 


38.4 


45.6 


16.5 


MOCASSIN model 


3.03 


4.1 


29.0 


57.0 


13.3 


H88 model 


0.95 


2.63 


40.0 


47.1 


11.0 



Table 6. Dust attenuation of resonance line radiation. The line flux intensi- 
ties, 1( A) are given relative to H/3 on a scale where I(H/3) = 100. 



Percent Destroyed 1(A) (attenuated) 

MOC. H88 MOC. H88 Obs 





Thick 


Thin 


Thick 


Thin 








Lya 1216 


68 


44 


80 


23 


859 


728 


* 


Civ 1550 


54 


38 


45 


21 


1025 


950 


512 


N v 1240 


39 


32 


25 


20 


24 


52 


46 


Mg II 2800 


48 


18 


32 




42 


29 


* 


Si IV 1400 


65 


30 


40 


13 


9.1 


5.7 


9 



4.3.1 The emergent spectral energy distribution 

The SEDs from our model were convolved with the IRAS broad- 
band filter transmission functions and compared with the line- and 
colour-corrected IRAS Point Source Catalog fluxes. A full descrip- 
tion of the line and colour corrections applied is given by H88. 
Following H88, the Cohen & Barlow (1980) flux at 9.6^m is to 
be preferred to the IRAS 12/im flux, given the large ionic emission 
line contribution to the latter. Table[5]lists MOCASSIN'S model pre- 
dictions in the second row and the observed fluxes in the first row, 
H88's model predictions (row 3) are also given for comparison. 

Although the IR flux distribution is slightly different from H88 
for the 5 wavelength regions, the total IR (~5-300/im) luminosities 
are similar for the two codes (within 5%). The remaining differ- 
ences are however acceptable in light of the different approach em- 
plyed by our code with regards to the dust model, where we treat 
the individual grain sizes independently in the RT, as well in the 
final temperature and SED calculations, and also with regards to 
the RT between the two optical depth regions that is treated self- 
consistently in our 3D simulation. 

4.3.2 Dust attenuation of UV resonance lines 

The emergent resonance emission line intensities are heavily atten- 
uated by dust absorption. The escape probability method adopted 
for the RT of the resonance lines has already been described in 
Section |2~TT1 Table |6] lists MOCASSIN'S predictions for the dust- 
attenuated resonance line intensities and percentage destruction 
fractions in each of the two sectors. The observed values and the 
line intensities and percentage destruction fractions predicted by 
H88 for the two sectors are also listed for comparison. All line in- 
tensities are given relative to H/3, on a scale where I(H/3) = 100. 

There are evident differences between our results and those of 
H88. These are not surprising and are due to a variety of factors. 
The atomic data set available to MOCASSIN has naturally been ex- 
panded in the 17 years since the publication of H88's results, The 
only meaningful comparison between the two codes can be made 
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Table 7. Percentage dust heating contributions. 



Energy source 




Percentage Heatin 


I Contribution 






Thick 


Thin 


Total 




MOC. 


H88 


MOC. 


H88 


MOC. 


H88 


Lya 1216 


55.7 


67.6 


24.2 


17.8 


47.1 


57.3 


UV continuum 


15.9 


17.2 


39.4 


58.5 


22.5 


25.8 


C IV 1550 


26.0 


15.0 


35.3 


21.8 


28.6 


16.4 


Other lines 


2.4 


0.2 


1.1 


1.9 


1.8 


0.5 



in the case of Lya, here we see a reasonable agreement (to within 
12%) in the predicted total attenuated flux for this line. Even for 
this line there are differences in the predicted destruction fractions, 
with MOCASSIN predicting a lower fraction in the thick sector and 
a higher fraction in the thin sector. This is easily understood by 
considering that in our 3D model resonance line packets emitted 
in the thin sector can cross over to the thicker one, with a higher 
probability of being absorbed by the dust in the thick sector, and 
vice-versa for lines emitted in the thick sector. 

The infrared energy balance is summarised in Table ?? which 
lists the major contributions to dust heating. In agreement with 
H88's conclusions we also found that Lya is the dominant source 
of heating in the optically thick region, followed by UV continuous 
radiation. In the thin sector, where the lower gas denisties naturally 
result in lower line emission and where the destruction factors are 
much reduced, Lya provides about 24% of the total grain heating, 
with the continuous radiation becoming the dominant source here. 
Heating by absorption of the other resonance line photons is also 
significant in both sectors. The small dicrepancies between our re- 
sults and H88's can once again be explained in terms of the 3D 
transfer as well as differences in the atomic data sets. 



4.3.3 The dust temperatures 

Grain temperatures for each grain size are determined from the 
equation of radiative equilibrium, as discussed in Section l2~3l Fig- 
ure [3] shows the radius-dependent temperature distributions for 
grains of radius 0.04/im (hottest), 0.1//m and 0.40/im (coolest). 
Each data point represents the temperature of a grain in a given 
cell and all cells belonging to the optically thin density phase are 
plotted in the top panel, whilst those belonging to the optically thick 
bicone are plotted in the bottom panel. The difference in the grain 
temperatures in cells belonging to the two different optical depth 
regions is evident, with the more opaque regions being naturally 
hotter. The difference is further accentuated by the extra heating 
provided by the resonance line photons, a process that is more ef- 
ficient in the optically thick sector than in the thin one. The scatter 
of the data is only partially due to numerical noise; it is also a re- 
sult of radiation transfer between the optically thick cones and the 
surrounding thinner regions. 

While the grain temperatures in the two sectors seem to be in 
agreement with those shown in figure 3 of H88, a closer inspection 
reveals that our temperatures in the thick sector are consistently 
lower, whilst those in the thin sector are consistently higher than 
predicted by H88. The very small discrepancies are due to the fact 
that in the case of H88's ID composite models the resonance line 
emission produced in one sector will only be able to heat the grains 
within the same sector. In our 3D model, which is self-consistent at 
the boundaries between the two sectors, packets can be transfered 
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Figure 3. Radial grain temperatures for grains of radius 0.04/im (hottest), 
OAfim (medium) and 0.40/xm (coolest). The cells in thin sectors are plotted 
in the top panel and those from the thick sector are in the bottom panel. This 
figure is comparable to figure 3 of H88. 



throughout the volume and, therefore, dust grains in the optically 
thin region are also heated by resonance line packets produced in 
the thick biconical regions. Naturally, the converse is also true, but 
the density imbalance between the two sectors practically results 
in dust-heating leaking from the optically thick to the surrounding 
optically thin region. All our model results are consistent with this 
interpretation. 



5 CONCLUSIONS 

We have presented a new tool for the analysis of dusty photoionised 
plasma. The dust-enhanced version of the fully 3D photoionisation 
code MOCASSIN is suitable for the interpretation of observational 
data from such environments, including emission line spectra, ther- 
mal SEDs and narrow or broad-band images at optical or IR wave- 
lengths. 

Our results demonstrate the reliability of the new version of 
MOCASSIN that employs a MC approach to solving the RT prob- 
lem for dust and gas simultaneously. The benchmark tests showed 
that MOCASSIN was able to reproduce DUSTY's solutions for grain 
temperatures and SEDs for a number of spherically symmetric 
pure-dust benchmark models (Ivezic et al., 1997). The code also 
achieved a good level of accuracy in matching the results of a num- 
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ber of well established 2D and 3D dust RT codes for a set of 2D 
disk benchmarks proposed by Pascucci et al. (2004). 

A fully self-consistent 3D photoionisation and dust model was 
constructed for the PN NGC 3918, which was treated as a realistic 
dust and gas benchmark case, thanks to the existence of a previous 
ID dust and gas photoionisation model by H88. Our results were 
shown to be in reasonable agreement with those obtained by H88, 
within the limits imposed by the fundamental differences between 
the two approaches and the latter' s need to adopt a composite model 
of two volume weighted ID models to describe a 3D structure. 

The new version of MOCASSIN (Version 2.0) is now fully 
functional and it will be made freely available to the astronomical 
community shortly after the publication of this article 2 . 
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